/*
 * 27.cpp
 *
 * Demonstrating three solutions to exercise 1.1.27 of
 * Algorithms, Fourth Edition.
 *
 */

#include <math.h>
#include <cstdlib>

#include <chrono>
#include <string>


namespace algs4
{


static unsigned long long calltimes = 0;


static long double binomial(int N, int K, long double p)
{
    ++calltimes;
    if (N == 0 && K == 0) return 1.0L;
    if (N < 0 || K < 0) return 0.0L;
    return (1.0L-p) * binomial(N-1, K, p) + p * binomial(N-1, K-1, p);
}


static long double binomial1_impl(long double** b, int N, int K, long double p)
{
    if (N == 0 && K == 0) return 1.0L;
    if (N < 0 || K < 0) return 0.0L;
    if (b[N][K] == -1) {
        ++calltimes;
        b[N][K] = (1.0L-p) * binomial1_impl(b, N-1, K, p) + p * binomial1_impl(b, N-1, K-1, p);
    }
    return b[N][K];
}


static long double binomial1(int N, int K, long double p)
{
    long double x;

    long double** b = new long double*[N+1];
    long double* data = new long double[(N+1)*(K+1)];

    int i, j;
    for (i = 0; i <= N; ++i)
        b[i] = data + (K+1) * i;
    for (i = 0; i <= N; ++i)
        for (j = 0; j <= K; ++j)
            b[i][j] = -1;

    x = binomial1_impl(b, N, K, p);

    delete[] data;
    delete[] b;

    return x;
}


static long double binomial2(int N, int K, long double p)
{
    long double x;

    if (N < K) {

        x = 0.0L;

    } else if (N == K) {

        x = powl(p, N);

    } else {

        int i, j;
        long double** b = new long double*[N+1];
        long double* data = new long double[(N+1)*(K+1)];
        for (i = 0; i <= N; ++i)
            b[i] = data + (K+1) * i;

        b[0][0] = 1.0L;

        // process '#' area
        // calcuate [1..N-K][0]
        for (i = 1; i <= N-K; ++i)
            b[i][0] = powl(1.0L-p, i);
        // calcuate [0][1..K]
        for (j = 1; j <= K; ++j)
            b[0][j] = 0.0L;
        // calcuate [1..N-K][1..K]
        for (i = 1; i <= N-K; ++i)
            for (j = 1; j <= K; ++j)
                b[i][j] = (1.0L-p) * b[i-1][j] + p * b[i-1][j-1];

        // process '*' area
        for (i = N-K+1; i <= N; ++i)
            for (j = i-(N-K); j <= K; ++j)
                b[i][j] = (1.0L-p) * b[i-1][j] + p * b[i-1][j-1];

        x = b[N][K];

        delete[] data;
        delete[] b;

    }

    return x;
}


static long double combination(int n, int m)
{
    if (m > n/2)
        m = n - m;
    long double a = 1L, b = 1L;
    int cnt = m, i;
    for (i = 1; i <= cnt; ++i) {
        b *= n--;
        a *= m--;
    }
    return b / a;
}


static long double binomial3(int N, int K, long double p)
{
    long double x;

    if (N < K) {

        x = 0.0L;

    } else if (N == K) {

        x = powl(p, N);

    } else {

        x = combination(N, K) * powl(1.0L-p, N-K) * powl(p, K);

    }

    return x;
}


} // end of namespace algs4


int main(int argc, char* argv[])
{
    using std::chrono::steady_clock;
    using std::chrono::time_point;
    using std::chrono::microseconds;

    time_point<steady_clock> start, end;

    if (argc != 4) {
        std::printf("\nUsage: /path/to/27.out <integer> <integer> <double>\n");
        std::exit(EXIT_FAILURE);
    }
    std::printf("\n");

    int N = std::stoi(argv[1]), K = std::stoi(argv[2]);
    long double p = std::stold(argv[3]), x;

    for (int i = 1; i <= 10; ++i) {
        std::printf("\n############################# Round %2d #\n", i);

        start = steady_clock::now();
        algs4::calltimes = 0;
        x = algs4::binomial(N, K, p);
        end = steady_clock::now();
        std::printf("binomial(%d, %d, %.2Lf):  %.12Lf\n"
                        "\tCalled %llu times.\n"
                        "\tTook %llu microseconds.\n\n",
                        N, K, p, x,
                        algs4::calltimes,
                        std::chrono::duration_cast<microseconds>(end - start).count());

        start = steady_clock::now();
        algs4::calltimes = 0;
        x = algs4::binomial1(N, K, p);
        end = steady_clock::now();
        std::printf("binomial1(%d, %d, %.2Lf): %.12Lf\n"
                        "\tCalled %llu times.\n"
                        "\tTook %llu microseconds.\n\n",
                        N, K, p, x,
                        algs4::calltimes,
                        std::chrono::duration_cast<microseconds>(end - start).count());

        start = steady_clock::now();
        x = algs4::binomial2(N, K, p);
        end = steady_clock::now();
        std::printf("binomial2(%d, %d, %.2Lf): %.12Lf\n"
                        "\tTook %llu microseconds.\n\n",
                        N, K, p, x,
                        std::chrono::duration_cast<microseconds>(end - start).count());

        start = steady_clock::now();
        x = algs4::binomial3(N, K, p);
        end = steady_clock::now();
        std::printf("binomial3(%d, %d, %.2Lf): %.12Lf\n"
                        "\tTook %llu microseconds.\n\n",
                        N, K, p, x,
                        std::chrono::duration_cast<microseconds>(end - start).count());

    }

    std::printf("\n");
    return 0;
}

